############################################################
##         Replication Code for Online Appendix E         ##
############################################################

###########################################################
##                 Figures OA.E1 & OA.E2                 ##
##                Roll Call Vote Samples                 ##
###########################################################

rm(list=ls())
library(superheat)
library(RColorBrewer)

#Figure OA.E1
alpha <- .5 
k <- 0.1 #0.2 #0.3

#Figure OA.E2
#alpha <- .25 #.5 #.75
#k <- 0.2

N <- 100
phi <- seq(.2,.8,.05)
nsims <- 50

#bills and sqs
b1 <- seq(0,1,.1)
b2 <- b1
sq1 <- seq(0,1,.1)
sq2 <- sq1
v <- length(b1)*length(sq1)

#heterogeneity scores
d1 <- seq(0.5,2, .1)
d2 <- d1

mat.percent <- matrix(NA,nrow=length(phi),ncol=length(d1))
mat.rcvcount <- matrix(NA,nrow=length(phi),ncol=length(d1))
mat.propcount <- matrix(NA,nrow=length(phi),ncol=length(d1))
percent.array <- array(NA,dim=c(length(phi),ncol=length(d1),nsims))

for(L in 1:(nsims)){
  for(K in 1:length(phi)){  
  for(J in 1:length(d1)){
      npL <- round((1-phi[K])*N,digits=0)
      npR <- round((phi[K])*N,digits=0)

      hypo.pro.1d <- cbind(rep(1,length(b1)*length(sq1)),rep(d1[J],length(b1)*length(sq1)),rep(phi[K],length(b1)*length(sq1)),
                           rep(b1,length(sq1)),rep(sq1,each=length(b1)))
      colnames(hypo.pro.1d) <- c("dim","d","phi","b","sq")
      hypo.pro.1d <- as.data.frame(hypo.pro.1d)
      hypo.pro.2d <- cbind(rep(2,length(b2)*length(sq2)),rep(d2[J],length(b2)*length(sq2)),rep(phi[K],length(b2)*length(sq2)),
                           rep(b2,length(sq2)),rep(sq2,each=length(b2)))
      colnames(hypo.pro.2d) <- c("dim","d","phi","b","sq")
      hypo.pro.2d <- as.data.frame(hypo.pro.2d)
      hypo.pro.all <- as.data.frame(rbind(hypo.pro.1d,hypo.pro.2d))
      hypo.pro.all$id <- seq(1,nrow(hypo.pro.all),1)
      
      #ideal points
      xL1 <- c(0,runif(npL-1,min=0-d1[J],max=d1[J])) #Going to use this that i've put the party leaders in the first and last spots
      xL2 <- xL1#c(0,runif(npL-1,min=0-d2[k],max=d2[k]))
      xR1 <- c(runif(npR-1,min=1-d1[J],max=1+d1[J]),1) #Going to use this that i've put the party leaders in the first and last spots
      xR2 <- xR1#c(runif(npR-1,min=1-d2[k],max=1+d2[k]),1)
      
      #matrices of legislators ideal points
      left.party <- as.data.frame(cbind(xL1,xL2,rep(0,npL)))
      colnames(left.party) <- c("x1","x2","partyid")
      right.party <- as.data.frame(cbind(xR1,xR2,rep(1,npR)))
      colnames(right.party) <- c("x1","x2","partyid")
      all.leg <- as.data.frame(rbind(left.party,right.party))
      
      #who votes for it?
      vote.func1 <- function(x1,b1,sq1){
        as.numeric(abs(x1-b1)  < abs(x1-sq1) )
      }
      vote.func2 <- function(x2,b2,sq2){
        as.numeric(abs(x2-b2)  < abs(x2-sq2) )
      }
      
      #full voting matrix (exogenous bills)
      out <- matrix(data=NA,nrow=N,ncol=v*2);dim(out)
      for(i in 1:N){
        for(j in 1:v){
          out[i,j] <- vote.func1(all.leg$x1[i],hypo.pro.1d$b[j],hypo.pro.1d$sq[j])
        }
      }
      for(i in 1:N){
        for(j in (v+1):(2*v)){
          out[i,j] <- vote.func2(all.leg$x2[i],hypo.pro.2d$b[j-v],hypo.pro.2d$sq[j-v])
        }
      }
      
      #does it pass?
      passes <- as.numeric(colSums(out) > .5*N)
      hypo.pro.all <- as.data.frame(cbind(hypo.pro.all,passes))
      for (i in 1:nrow(hypo.pro.all)){
        hypo.pro.all$outcome[i] <- if(hypo.pro.all$passes[i]==0) hypo.pro.all$sq[i] else hypo.pro.all$b[i]
      }
      
      #cohesion scores by legislator
      left.aligned <- rep(NA,length=ncol(out))
      for(i in 1:ncol(out)){
        left.aligned[i] <- round((sum(as.numeric(out[2:npL,i]==out[1,i])))/npL,4)
      }
      right.aligned <- rep(NA,length=ncol(out))
      for(i in 1:ncol(out)){
        right.aligned[i] <- round((sum(as.numeric(out[(npL+1):(N-1),i]==out[N,i])))/npR,4)
      }
      out.cohesion <- matrix(data=NA,nrow=N,ncol=v*2)
      for(j in 1:(v*2)){
        for(i in 1:npL){
          out.cohesion[i,j] <- if (out[i,j] == out[1,j]) {
            left.aligned[j] - right.aligned[j]
          } else {
            0
          }
        }
      }
      for(j in 1:(v*2)){
        for(i in (npL+1):N){
          out.cohesion[i,j] <- if (out[i,j] == out[N,j]) {
            right.aligned[j]-left.aligned[j]
          } else {
            0
          }
        }
      }
      out.cohesion.abs <- matrix(data=NA,nrow=N,ncol=v*2)
      for(j in 1:(v*2)){
        for(i in 1:npL){
          out.cohesion.abs[i,j] <- if (out[i,j] == out[1,j]) {
            left.aligned[j] 
          } else {
            0
          }
        }
      }
      for(j in 1:(v*2)){
        for(i in (npL+1):N){
          out.cohesion.abs[i,j] <- if (out[i,j] == out[N,j]) {
            right.aligned[j]
          } else {
            0
          }
        }
      }

    #roll call vote requests?
    out.cohesion.kmat <- matrix(data=NA,nrow=N,ncol=v*2)
      for(i in 1:N){
        for(j in 1:(v*2)){
          out.cohesion.kmat[i,j] <- as.numeric(out.cohesion[i,j] > k)
        }
      }
      hypo.pro.all$left.cscore <- out.cohesion[1,] 
      hypo.pro.all$right.cscore <- out.cohesion[N,]
      hypo.pro.all$left.cscore.abs <- out.cohesion.abs[1,] 
      hypo.pro.all$right.cscore.abs <- out.cohesion.abs[N,]
      hypo.pro.all$left.rcv.request <- out.cohesion.kmat[1,]
      hypo.pro.all$right.rcv.request <- out.cohesion.kmat[N,]
      hypo.pro.all$rcv.request <- as.numeric(colSums(out.cohesion.kmat)>0)
      
    #proposed?
      for(i in 1:nrow(hypo.pro.all)){
        hypo.pro.all$proposing.party[i] <- if(hypo.pro.all$b[i] < hypo.pro.all$sq[i]) 0 else 1 
      }
      hypo.pro.all$proposed.left <- rep(0,length(nrow(hypo.pro.all)))
      for(i in 1:nrow(hypo.pro.all)){
        if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
          hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$b[i])+hypo.pro.all$left.cscore[i]-k > alpha*-abs(0-hypo.pro.all$sq[i]))
        } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
          hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$sq[i])+hypo.pro.all$left.cscore[i]-k > alpha*-abs(0-hypo.pro.all$sq[i]))
        } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==1) {
          hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$b[i]) > alpha*-abs(0-hypo.pro.all$sq[i]))
        } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==0) {
          hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$sq[i]) > alpha*-abs(0-hypo.pro.all$sq[i]))
        } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
          hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$b[i])+hypo.pro.all$left.cscore[i] > alpha*-abs(0-hypo.pro.all$sq[i]))
        } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
          hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$sq[i])+hypo.pro.all$left.cscore[i] > alpha*-abs(0-hypo.pro.all$sq[i]))
        } else {
          0
        }
      }
      hypo.pro.all$proposed.right <- rep(0,length(nrow(hypo.pro.all)))
      for(i in 1:nrow(hypo.pro.all)){
        if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
          hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$b[i])+hypo.pro.all$right.cscore[i]-k > alpha*-abs(1-hypo.pro.all$sq[i]))
        } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
          hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$sq[i])+hypo.pro.all$right.cscore[i]-k > alpha*-abs(1-hypo.pro.all$sq[i]))
        } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==1) {
          hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$b[i]) > alpha*-abs(1-hypo.pro.all$sq[i]))
        } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==0) {
          hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$sq[i]) > alpha*-abs(1-hypo.pro.all$sq[i]))
        } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
          hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$b[i])+hypo.pro.all$right.cscore[i] > alpha*-abs(1-hypo.pro.all$sq[i]))
        } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
          hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$sq[i])+hypo.pro.all$right.cscore[i] > alpha*-abs(1-hypo.pro.all$sq[i]))
        } else {
          0
        }
      }
      hypo.pro.all$proposed <- hypo.pro.all$proposed.left+hypo.pro.all$proposed.right
      all.proposals <- hypo.pro.all[which(hypo.pro.all$proposed==1),]
      rcv.proposals <- hypo.pro.all[which(hypo.pro.all$proposed==1 & hypo.pro.all$rcv.request==1),]

    mat.percent[K,J] <- nrow(rcv.proposals) / nrow(all.proposals)
    mat.rcvcount[K,J] <- nrow(rcv.proposals)
    mat.propcount[K,J] <- nrow(all.proposals)
    percent.array[K,J,L] <- mat.percent[K,J]
    }
  }
}

###########################################################
#content for labeling heatmap figures

  rcv.plot <- apply(percent.array,1:2,mean)
  d.lab <- c("0.5","","","","","1.0","","","","","1.5","","","","","2")
  heat.cols <- RColorBrewer::brewer.pal(10, "Greys")
  rownames(rcv.plot) <- c("","0.25","","","","","0.5","","","","","0.75","")#phi
  colnames(rcv.plot) <- c("0.5","","","","","1","","","","","1.5","","","","","2")

#pdf(".pdf",width=8,height=9)
superheat(X=rcv.plot,heat.lim=c(0,1),
          heat.pal=heat.cols,heat.na.col = "white",
          grid.hline.col = "white",grid.vline.col = "white",
          title="alpha = 0.X, k = 0.X", title.size=6,
          column.title = "Party Heterogeneity\n (d)", row.title.size = 6,
          row.title = "Right Party Seat Share\n (phi)  ", column.title.size = 6,
          left.label.size = .1,left.label.col = "white", left.label.text.size = 5,
          bottom.label.size = .1,bottom.label.col = "white", bottom.label.text.size = 5)
#dev.off()

###########################################################
##                      Table OA.E1                      ##
###########################################################

#setwd("/Main Text-Figures/Figure2/")
load("idealpoints-fig2-data.RData")

#For Full Sample
  #row 1
  summary(as.numeric(table(full.min.plot)))
  sd(as.numeric(table(full.min.plot)))
  #row 2
  summary(as.numeric(table(full.med.plot)))
  sd(as.numeric(table(full.med.plot)))
  #row 3
  summary(as.numeric(table(full.max.plot)))
  sd(as.numeric(table(full.max.plot)))

#For RCV Sample
  #row 4
  summary(as.numeric(table(rcv.min.plot)))
  sd(as.numeric(table(rcv.min.plot)))
  #row 5
  summary(as.numeric(table(rcv.med.plot)))
  sd(as.numeric(table(rcv.med.plot)))
  #row 6
  summary(as.numeric(table(rcv.max.plot)))
  sd(as.numeric(table(rcv.max.plot)))

###########################################################
##                     Figure OA.E3                      ##
###########################################################

#see code from "Main Text-Figures/Figure1/idealpoints-fig1-code.R"

###########################################################
##                     Figure OA.E4                      ##
###########################################################

#see code from "Main Text-Figures/Figure2/idealpoints-fig2-code.R"

###########################################################
##                     Figure OA.E5                      ##
###########################################################

#see code from "Main Text-Figures/Figure1/idealpoints-fig1-code.R"

###########################################################
##                     Figure OA.E6                      ##
###########################################################

#see code from "Main Text-Figures/Figure1/idealpoints-fig1-code.R"

###########################################################
##                     Figure OA.E7                      ##
###########################################################

alpha <- .5
k <- .25
N <- 100
v <- 500
nsim <- 15
nsims <- 1

prev.all <- NULL
rcv.count.all <- NULL
rcv.median.all <- NULL
full.median.all <- NULL

for(K in 1:nsims){
prev <- NULL
rcv.count <- NULL
rcv.median <- NULL
full.median <- NULL

for(J in 1:nsim){
  phi <- runif(1,.2,.8)
  b1 <- runif(v,0,1)
  b2 <- b1#runif(v,0,1)
  sq1 <- runif(v,0,1)
  sq2 <- sq1#runif(v,0,1)
  npL <- round((1-phi)*N,digits=0)
  npR <- round((phi)*N,digits=0)
  
  d1 <- runif(1,.5,2)
  d2 <- d1#runif(1,.5,2)

  hypo.pro.1d <- cbind(rep(1,length(b1)),rep(d1,length(b1)),b1,sq1)
  colnames(hypo.pro.1d) <- c("dim","d","b","sq")
  hypo.pro.2d <- cbind(rep(2,length(b2)),rep(d2,length(b2)),b2,sq2)
  colnames(hypo.pro.2d) <- c("dim","d","b","sq")
  hypo.pro.all <- as.data.frame(rbind(hypo.pro.1d,hypo.pro.2d))
  hypo.pro.all$id <- seq(1,nrow(hypo.pro.all),1)
  
  ########################################

  xL1 <- c(0,runif(npL-1,min=0-d1,max=d1)) #Going to use this that i've put the party leaders in the first and last spots
  xL2 <- xL1#c(0,runif(npL-1,min=0-d2,max=d2))
  xR1 <- c(runif(npR-1,min=1-d1,max=1+d1),1) #Going to use this that i've put the party leaders in the first and last spots
  xR2 <- xR1#c(runif(npR-1,min=1-d2,max=1+d2),1)
  left.party <- as.data.frame(cbind(round(xL1,3),round(xL2,3),rep(0,npL)))
  colnames(left.party) <- c("x1","x2","partyid")
  right.party <- as.data.frame(cbind(round(xR1,3),round(xR2,3),rep(1,npR)))
  colnames(right.party) <- c("x1","x2","partyid")
  all.leg <- as.data.frame(rbind(left.party,right.party))
  
  ########################################
  #who votes for it?
  vote.func1 <- function(x1,b1,sq1){
    as.numeric(abs(x1-b1)  < abs(x1-sq1) )
  }
  vote.func2 <- function(x2,b2,sq2){
    as.numeric(abs(x2-b2)  < abs(x2-sq2) )
  }
  out <- matrix(data=NA,nrow=N,ncol=v*2);dim(out)
  for(i in 1:N){
    for(j in 1:v){
      out[i,j] <- vote.func1(all.leg$x1[i],b1[j],sq1[j])
    }
  }
  for(i in 1:N){
    for(j in (v+1):(2*v)){
      out[i,j] <- vote.func2(all.leg$x2[i],b2[j-v],sq2[j-v])
    }
  }
  
  #passes?
  passes <- as.numeric(colSums(out) > .5*N)
  hypo.pro.all <- as.data.frame(cbind(hypo.pro.all,passes))
  for (i in 1:nrow(hypo.pro.all)){
    hypo.pro.all$outcome[i] <- if(hypo.pro.all$passes[i]==0) hypo.pro.all$sq[i] else hypo.pro.all$b[i]
  }
  left.aligned <- rep(NA,length=ncol(out))
  for(i in 1:ncol(out)){
    left.aligned[i] <- round((sum(as.numeric(out[2:npL,i]==out[1,i])))/npL,4)
  }
  right.aligned <- rep(NA,length=ncol(out))
  for(i in 1:ncol(out)){
    right.aligned[i] <- round((sum(as.numeric(out[(npL+1):(N-1),i]==out[N,i])))/npR,4)
  }
  out.cohesion <- matrix(data=NA,nrow=N,ncol=v*2)
  for(j in 1:(v*2)){
    for(i in 1:npL){
      out.cohesion[i,j] <- if (out[i,j] == out[1,j]) {
        left.aligned[j] - right.aligned[j]
      } else {
        0
      }
    }
  }
  for(j in 1:(v*2)){
    for(i in (npL+1):N){
      out.cohesion[i,j] <- if (out[i,j] == out[N,j]) {
        right.aligned[j]-left.aligned[j]
      } else {
        0
      }
    }
  }
  out.cohesion.abs <- matrix(data=NA,nrow=N,ncol=v*2)
  for(j in 1:(v*2)){
    for(i in 1:npL){
      out.cohesion.abs[i,j] <- if (out[i,j] == out[1,j]) {
        left.aligned[j] 
      } else {
        0
      }
    }
  }
  for(j in 1:(v*2)){
    for(i in (npL+1):N){
      out.cohesion.abs[i,j] <- if (out[i,j] == out[N,j]) {
        right.aligned[j]
      } else {
        0
      }
    }
  }
  ########################################
  #roll call votes requested?
  out.cohesion.kmat <- matrix(data=NA,nrow=N,ncol=v*2)
  for(i in 1:N){
    for(j in 1:(v*2)){
      out.cohesion.kmat[i,j] <- as.numeric(out.cohesion[i,j] > k)
    }
  }
  hypo.pro.all$left.cscore <- out.cohesion[1,] 
  hypo.pro.all$right.cscore <- out.cohesion[N,]
  hypo.pro.all$left.cscore.abs <- out.cohesion.abs[1,] 
  hypo.pro.all$right.cscore.abs <- out.cohesion.abs[N,]
  hypo.pro.all$left.rcv.request <- out.cohesion.kmat[1,]
  hypo.pro.all$right.rcv.request <- out.cohesion.kmat[N,]
  hypo.pro.all$rcv.request <- as.numeric(colSums(out.cohesion.kmat)>0)
  all.proposals <- hypo.pro.all
  rcv.proposals <- hypo.pro.all[which(hypo.pro.all$rcv.request==1),]
  
  ########################################
  #vote matrixes of proposed votes
  vote.matrix.all <- out[,all.proposals$id]
  vote.matrix.rcv <- out[,rcv.proposals$id]
  
  ########################################
  #ideal point estimation
  y.all <- as.data.frame(vote.matrix.all)
  rc.all <- rollcall(y.all)
  rc.all$votes[rc.all$votes==0]<- -1
  rc.all$codes$nay <- -1
  rcfast.all <- convertRC(rc.all)
  p.all <- makePriors(rcfast.all$n, rcfast.all$m, 1)
  s.all <- getStarts(rcfast.all$n, rcfast.all$m, 1)
  fast.out.all <- binIRT(.rc = rcfast.all,
                         .starts = s.all,
                         .priors = p.all,
                         .control = {
                           list(threads = 6,
                                verbose = FALSE,
                                thresh = 1e-6
                           )
                         },
                         .anchor_subject=100
  )
  full.est.all <- round(fast.out.all$means$x,3)
  full.est.all.party <- cbind(seq(1:length(fast.out.all$means$x)),full.est.all,all.leg)
  colnames(full.est.all.party) <- c("id","full.est.irt","idp1.true","idp2.true","party.id")
  
  y.rcv <- as.data.frame(vote.matrix.rcv)
  rc.rcv <- rollcall(y.rcv)
  rc.rcv$votes[rc.rcv$votes==0]<- -1
  rc.rcv$codes$nay <- -1
  rcfast.rcv <- convertRC(rc.rcv)
  p.rcv <- makePriors(rcfast.rcv$n, rcfast.rcv$m, 1)
  s.rcv <- getStarts(rcfast.rcv$n, rcfast.rcv$m, 1)
  fast.out.rcv <- binIRT(.rc = rcfast.rcv,
                         .starts = s.rcv,
                         .priors = p.rcv,
                         .control = {
                           list(threads = 6,
                                verbose = FALSE,
                                thresh = 1e-6
                           )
                         },
                         .anchor_subject=100
  )
  full.est.rcv <- round(fast.out.rcv$means$x,3)
  full.est.rcv.party <- cbind(seq(1:length(fast.out.rcv$means$x)),full.est.rcv,all.leg)
  colnames(full.est.rcv.party) <- c("id","rcv.est.irt","idp1.true","idp2.true","party.id")

  ########################################
  dat <- full.est.rcv.party
  dat$full.est.irt <- full.est.all.party$full.est.irt
  colnames(dat) <- c("id","partial","true","true2","party","full")
  dat.order <- dat[order(dat$true),]
  dat.order <- as.data.frame(cbind(seq(1:N),dat.order))
  colnames(dat.order) <- c("true.rank","id","partial","true","true2","party","full")
  dat.left <- dat.order[which(dat.order$party==0),]
  dat.right <- dat.order[which(dat.order$party==1),]
  
  s.prev <- ncol(y.rcv)/ncol(y.all)
  s.rcv <- ncol(y.rcv)
  s.rcv.median <- max(100-as.numeric(table(dat$partial == median(dat$partial,na.rm=TRUE)))[1],1)
  s.full.median <- max(100-as.numeric(table(dat$full == median(dat$full,na.rm=TRUE)))[1],1)
  prev <- c(prev,s.prev)
  rcv.count <- c(rcv.count,s.rcv)
  rcv.median <- c(rcv.median,s.rcv.median)
  full.median <- c(full.median,s.full.median)
}  

prev.all <- c(prev.all,mean(prev))
rcv.count.all <- c(rcv.count.all, mean(rcv.count))
rcv.median.all <- c(rcv.median.all,mean(rcv.median))
full.median.all <- c(full.median.all,mean(full.median))
}

ex.dat <- as.data.frame(cbind(prev.all,rcv.count.all,rcv.median.all,full.median.all))
ex.dat <- ex.dat[order(ex.dat$prev.all),]

###########################################################
#Generate plot OA.E7
#setwd("/OA.E/")
#load("fig-oae7-data.RData")

#pdf(".pdf",width=6,height=5)
  par(mar=c(4,4,1,1))
  plot(x=ex.dat$prev.all,y=ex.dat$rcv.median.all,pch=16,xlab="",ylab="",
     ylim=c(0,50),xlim=c(.1,.7),xaxt='n',yaxt='n')
  lines(x=ex.dat$prev.all,y=ex.dat$rcv.median.all,lwd=2)
  lines(x=ex.dat$prev.all,y=ex.dat$full.median.all,col="grey60")
  points(x=ex.dat$prev.all,y=ex.dat$full.median.all,pch=15,col="grey60")
  mtext(side=1,line=2,"Prevalence of Roll Call Votes")
  mtext(side=2,line=3,"% Legislators Co-Located at Chamber Median")
  axis(side=2,at=c(0,10,20,30,40,50),labels=c("0%","10%","20%","30%","40%","50%"),las=1)
  axis(side=1,at=c(.1,.2,.3,.4,.5,.6,.7,.8),labels=c("","20%","","40%","","60%","","80%"))
#dev.off()

###########################################################
##                      Table OA.E2                      ##
###########################################################
#setwd("/Main Text-Figures/Figure4/")
load("cohesion-fig4-data.RData")

#For Signed Difference
  #row 1
  summary(all);sd(all)
  #row 2
  summary(sim.info1$diff);sd(sim.info1$diff)
#For Absolute Difference 
  #row 3
  summary(abs(all));sd(abs(all))
  #row4
  summary(sim.info1$abs.diff);sd(sim.info1$abs.diff)

###########################################################
##                 Figure OA.E8 & OA.310                 ##
###########################################################

#see code from "Main Text-Figures/Figure3/cohesion-fig3ab-code.R"
#see code from "Main Text-Figures/Figure3/cohesion-fig3cd-code.R"

###########################################################
##                 Figure OA.E9 & OA.E11                 ##
###########################################################

#see code from "Main Text-Figures/Figure4/cohesion-fig4-code.R"

###########################################################
##                     Figure OA.E12                     ##
###########################################################

alpha <- 1
k <- 0.1
N <- 100
v <- 500
nleg <- 50
nsim <- 50

diff.left.cohesion.sim <- NULL
diff.right.cohesion.sim <- NULL
cohesion.right.all.sim <- NULL
cohesion.right.rcv.sim <- NULL
cohesion.left.all.sim <- NULL
cohesion.left.rcv.sim <- NULL
sim.info.all <- NULL

for(z in 1:nleg){
  phi <- runif(1,.5,.8)
  npL <- round((1-phi)*N,digits=0)
  npR <- round((phi)*N,digits=0)
  d1 <- runif(n=1,min=0.5,max=2)
  d2 <- d1#runif(n=1,min=0.5,max=2)

  #ideal points
  xL1 <- c(0,runif(npL-1,min=0-d1,max=d1)) #Going to use this that i've put the party leaders in the first and last spots
  xL2 <- xL1#c(0,runif(npL-1,min=0-d2[k],max=d2[k]))
  xR1 <- c(runif(npR-1,min=1-d1,max=1+d1),1) #Going to use this that i've put the party leaders in the first and last spots
  xR2 <- xR1#c(runif(npR-1,min=1-d2[k],max=1+d2[k]),1)
  left.party <- as.data.frame(cbind(xL1,xL2,rep(0,npL)))
  colnames(left.party) <- c("x1","x2","partyid")
  right.party <- as.data.frame(cbind(xR1,xR2,rep(1,npR)))
  colnames(right.party) <- c("x1","x2","partyid")
  all.leg <- as.data.frame(rbind(left.party,right.party))

  for(j in 1:nsim){
    b1 <- runif(v,0,1)
    b2 <- b1#runif(v,0,1)
    sq1 <- runif(v,0,1)
    sq2 <- sq1
  
    hypo.pro.1d <- cbind(rep(1,length(b1)),rep(d1,length(b1)),b1,sq1)
    colnames(hypo.pro.1d) <- c("dim","d","b","sq")
    hypo.pro.2d <- cbind(rep(2,length(b2)),rep(d2,length(b2)),b2,sq2)
    colnames(hypo.pro.2d) <- c("dim","d","b","sq")
    hypo.pro.all <- as.data.frame(rbind(hypo.pro.1d,hypo.pro.2d))
    hypo.pro.all$id <- seq(1,nrow(hypo.pro.all),1)
    
    #who votes for it?
    vote.func1 <- function(x1,b1,sq1){
      as.numeric(abs(x1-b1)  < abs(x1-sq1) )
    }
    vote.func2 <- function(x2,b2,sq2){
      as.numeric(abs(x2-b2)  < abs(x2-sq2) )
    }
    out <- matrix(data=NA,nrow=N,ncol=v*2);dim(out)
    
    for(i in 1:N){
      for(j in 1:v){
        out[i,j] <- vote.func1(all.leg$x1[i],b1[j],sq1[j])
      }
    }
    
    for(i in 1:N){
      for(j in (v+1):(2*v)){
        out[i,j] <- vote.func2(all.leg$x2[i],b2[j-v],sq2[j-v])
      }
    }
    
    #passes?
    passes <- as.numeric(colSums(out) > .5*N)
    hypo.pro.all <- as.data.frame(cbind(hypo.pro.all,passes))
    for (i in 1:nrow(hypo.pro.all)){
      hypo.pro.all$outcome[i] <- if(hypo.pro.all$passes[i]==0) hypo.pro.all$sq[i] else hypo.pro.all$b[i]
    }

    #cohesion scores?
    left.aligned <- rep(NA,length=ncol(out))
    for(i in 1:ncol(out)){
      left.aligned[i] <- round((sum(as.numeric(out[2:npL,i]==out[1,i])))/npL,4)
    }
    right.aligned <- rep(NA,length=ncol(out))
    for(i in 1:ncol(out)){
      right.aligned[i] <- round((sum(as.numeric(out[(npL+1):(N-1),i]==out[N,i])))/npR,4)
    }
    out.cohesion <- matrix(data=NA,nrow=N,ncol=v*2)
    for(j in 1:(v*2)){
      for(i in 1:npL){
        out.cohesion[i,j] <- if (out[i,j] == out[1,j]) {
          left.aligned[j] - right.aligned[j]
        } else {
          0
        }
      }
    }
    for(j in 1:(v*2)){
      for(i in (npL+1):N){
        out.cohesion[i,j] <- if (out[i,j] == out[N,j]) {
          right.aligned[j]-left.aligned[j]
        } else {
          0
        }
      }
    }
    out.cohesion.abs <- matrix(data=NA,nrow=N,ncol=v*2)
    for(j in 1:(v*2)){
      for(i in 1:npL){
        out.cohesion.abs[i,j] <- if (out[i,j] == out[1,j]) {
          left.aligned[j] 
        } else {
          0
        }
      }
    }
    for(j in 1:(v*2)){
      for(i in (npL+1):N){
        out.cohesion.abs[i,j] <- if (out[i,j] == out[N,j]) {
          right.aligned[j]
        } else {
          0
        }
      }
    }
    #roll call votes requested?
    out.cohesion.kmat <- matrix(data=NA,nrow=N,ncol=v*2)
    for(i in 1:N){
      for(j in 1:(v*2)){
        out.cohesion.kmat[i,j] <- as.numeric(out.cohesion[i,j] > k)
      }
    }
    hypo.pro.all$left.cscore <- out.cohesion[1,] 
    hypo.pro.all$right.cscore <- out.cohesion[N,]
    hypo.pro.all$left.cscore.abs <- out.cohesion.abs[1,] 
    hypo.pro.all$right.cscore.abs <- out.cohesion.abs[N,]
    hypo.pro.all$left.rcv.request <- out.cohesion.kmat[1,]
    hypo.pro.all$right.rcv.request <- out.cohesion.kmat[N,]
    hypo.pro.all$rcv.request <- as.numeric(colSums(out.cohesion.kmat)>0)
    all.proposals <- hypo.pro.all
    
    all.proposals <- all.proposals[which(all.proposals$rcv.request==0),]
    rcv.proposals <- hypo.pro.all[which(hypo.pro.all$rcv.request==1),]
    
    #vote matrixes of proposed votes
    vote.matrix.all <- out[,all.proposals$id]
    vote.matrix.rcv <- out[,rcv.proposals$id]        
    diff.left.cohesion <- mean(all.proposals$left.cscore.abs,na.rm=TRUE) - mean(rcv.proposals$left.cscore.abs,na.rm=TRUE)
    diff.right.cohesion <- mean(all.proposals$right.cscore.abs,na.rm=TRUE) - mean(rcv.proposals$right.cscore.abs,na.rm=TRUE)
  }

  diff.left.cohesion.sim <-  c(diff.left.cohesion.sim,diff.left.cohesion)
  diff.right.cohesion.sim <- c(diff.right.cohesion.sim,diff.right.cohesion)
  sim.info <- c(nrow(all.proposals),nrow(rcv.proposals),phi,d1)
  sim.info.all <- rbind(sim.info.all,sim.info)
}

sim.info.sum <- sim.info.all
abs.left <- abs(diff.left.cohesion.sim)
abs.right <- abs(diff.right.cohesion.sim)
abs.all <- c(abs.left,abs.right)
all <- c(diff.left.cohesion.sim,diff.right.cohesion.sim)

sim.info.left <- cbind(sim.info.all,diff.left.cohesion.sim)
colnames(sim.info.left) <- c("non.rcv.count","rcv.count","phi","d","diff")
sim.info.right <- cbind(sim.info.all,diff.right.cohesion.sim)
colnames(sim.info.right) <- c("non.rcv.count","rcv.count","phi","d","diff")
sim.info.all <- rbind(sim.info.left,sim.info.right)
sim.info.all <- as.data.frame(sim.info.all)
sim.info.all$abs.diff <- abs(sim.info.all$diff)

sim.info1 <- sim.info.all[which(sim.info.all$phi >=.35 & sim.info.all$phi <= .65 & sim.info.all$d >=.75 & sim.info.all$d <= 1.5),]
sim.info2 <- sim.info.all[which( (sim.info.all$phi< .35 | sim.info.all$phi > .65) & (sim.info.all$d < .75 | sim.info.all$d > 1.5)),]
sd(sim.info2$abs.diff)

###########################################################
#Generate plot OA.E12
#setwd("/OA.E/")
#load("fig-oae12-data.RData")


#pdf(".pdf",width=6,height=6)
par(mfrow=c(1,1),mar=c(5,5,3,2))
hist(all,border="white",col="grey60",
     xlab="Non-RCV Cohesion - RCV Cohesion",
     main="",xaxt='n',
     xlim=c(-.25,.25),ylim=c(0,12),yaxt='n',
     breaks=seq(-.5,.5,.025),ylab="")
#hist(sim.info1$diff,breaks=seq(-.5,.5,.025),border="white",col="grey60",add=T)
#abline(v=c(-.1,.1),lwd=2)
mtext(side=3,line=1,"k = 0.1",cex=1.5)
#segments(x0=.1,y0=0,x1=.1,y1=10,lwd=2)
#segments(x0=-.1,y0=0,x1=-.1,y1=10,lwd=2)
axis(side=2,at=seq(0,12,2),labels=seq(0,12,2),las=1)
axis(side=1,at=seq(-.25,.25,.05),labels=c("",-.2,"",-.1,"",0,"",.1,"",.2,""),las=1)
mtext(side=2,line=3,"Frequency")
#dev.off()

###########################################################
##                      Table OA.E3                      ##
##                          and                          ##
##                Figures OA.E13 & OA.E14                ##
###########################################################
#setwd("/OA.E/")

#column 1
#load("table-oae3-data-col1.RData")
mean(sim.info.all$total)
mean(sim.info.all$rcv)/mean(sim.info.all$total)
mean(sim.info.all2$percent.all.d1,na.rm=TRUE)
mean(sim.info.all2$percent.all.d2,na.rm=TRUE)
as.numeric(table(eigenall.order < 2.1))[1]*4
as.numeric(table(eigen.all[2,]<1))[1]*4 
mean(sim.info.all2$percent.rcv.d1,na.rm=TRUE)
mean(sim.info.all2$percent.rcv.d2,na.rm=TRUE)
as.numeric(table(eigenrcv.order < 2.1))[1]*4
as.numeric(table(eigen.rcv[2,]<1))[1]*4 

#Generate Figure OA.E13, Row 1
#pdf(".pdf",width=8,height=4)
par(mfrow=c(1,2),mar=c(4,4,4,1))
plot(eigen.all[,1],xlim=c(1,10),ylim=c(0,20),type='l',ylab="",xlab="",yaxt='n',col="grey60",lwd=1,xaxt='n')
abline(v=2,col="black",lwd=4)
for(i in 1:ncol(eigen.all)){
  lines(eigen.all[,i],col="grey60",lwd=1)
}
axis(side=2,at=seq(0,20,5),labels=seq(0,20,5),las=1)
mtext(side=3,line=2.25,"(a) All Votes",cex=1.5)
mtext(side=3,line=.75,print(paste("Average Total Votes  = ",round(mean(sim.info.all$total),0))),cex=1)
text(x=10,y=19,print(paste("k = ",k)),cex=1,pos=2)
text(x=10,y=15,print(paste("Missed Second\n Dimension: ",
                           100-as.numeric(table(eigen.all[2,]<1))[1]*4,"-",
                           100- as.numeric(table(eigenall.order < 2.1))[1]*4,"%")),pos=2)
mtext(side=2,line=3,"Eigenvalues")
mtext(side=1,line=2.5,"Number of Dimensions")
axis(side=1,at=seq(2,10,2),labels=seq(2,10,2))
axis(side=1,at=c(1,3),labels=c(1,3))
abline(h=1,lty=2,col="red")

par(mar=c(4,1,4,4))
plot(eigen.rcv[,1],xlim=c(1,10),pch=16,ylim=c(0,20),type='l',ylab="",xlab="",yaxt='n',col="grey60",lwd=1,xaxt='n')
abline(v=2,col="black",lwd=4)
for(i in 1:ncol(eigen.rcv)){
  lines(eigen.rcv[,i],col="grey60",lwd=1)
}
axis(side=4,at=seq(0,20,5),labels=seq(0,20,5),las=1)
mtext(side=3,line=2.25,"(b) Observed Roll Call Votes",cex=1.5)
mtext(side=3,line=.75,print(paste("Average % RCV  = ",round(mean(sim.info.all$percent.rcv),0),"%")),cex=1)
text(x=10,y=19,print(paste("k = ",k)),cex=1,pos=2)
text(x=10,y=15,print(paste0("Missed Second\n Dimension: ",
                            100-as.numeric(table(eigen.rcv[2,]<1))[1]*4,"-",
                            100-as.numeric(table(eigenrcv.order < 2.1))[1]*4,"%")),pos=2)
mtext(side=4,line=3,"Eigenvalues")
mtext(side=1,line=2.5,"Number of Dimensions")
axis(side=1,at=seq(2,10,2),labels=seq(2,10,2))
axis(side=1,at=c(1,3),labels=c(1,3))
abline(h=1,lty=2,col="red")
#dev.off()

#column 2
#load("Main Text-Figures/Figure5/dimensionality-fig5-data.RData")
sim.info.all <- sim.info.all[1:25,]
sim.info.all2 <- sim.info.all2[1:25,]
mean(sim.info.all$total)
mean(sim.info.all$rcv)/mean(sim.info.all$total)
mean(sim.info.all2$percent.all.d1,na.rm=TRUE)
mean(sim.info.all2$percent.all.d2,na.rm=TRUE)
as.numeric(table(eigenall.order[1:25] < 2))[1]*4
as.numeric(table(eigen.all[2,1:25]<1))[1]*4 
mean(sim.info.all2$percent.rcv.d1,na.rm=TRUE)
mean(sim.info.all2$percent.rcv.d2,na.rm=TRUE)
as.numeric(table(eigenrcv.order[1:25] < 2))[1]*4
as.numeric(table(eigen.rcv[2,1:25]<1))[1]*4 


#column 3
#load("table-oae3-data-col3.RData")
mean(sim.info.all$total)
mean(sim.info.all$rcv)/mean(sim.info.all$total)
mean(sim.info.all2$percent.all.d1,na.rm=TRUE)
mean(sim.info.all2$percent.all.d2,na.rm=TRUE)
as.numeric(table(eigenall.order < 2.1))[1]*4
as.numeric(table(eigen.all[2,]<1))[1]*4 
mean(sim.info.all2$percent.rcv.d1,na.rm=TRUE)
mean(sim.info.all2$percent.rcv.d2,na.rm=TRUE)
as.numeric(table(eigenrcv.order < 2.1))[1]*4
as.numeric(table(eigen.rcv[2,]<1))[1]*4 

#Generate Figure OA.E13, Row 2
#pdf(".pdf",width=8,height=4)
par(mfrow=c(1,2),mar=c(4,4,4,1))
plot(eigen.all[,1],xlim=c(1,10),ylim=c(0,20),type='l',ylab="",xlab="",yaxt='n',col="grey60",lwd=1,xaxt='n')
abline(v=2,col="black",lwd=4)
for(i in 1:ncol(eigen.all)){
  lines(eigen.all[,i],col="grey60",lwd=1)
}
axis(side=2,at=seq(0,20,5),labels=seq(0,20,5),las=1)
mtext(side=3,line=2.25,"(a) All Votes",cex=1.5)
mtext(side=3,line=.75,print(paste("Average Total Votes  = ",round(mean(sim.info.all$total),0))),cex=1)
text(x=10,y=19,print(paste("k = ",k)),cex=1,pos=2)
text(x=10,y=15,print(paste("Missed Second\n Dimension: ",
                           100-as.numeric(table(eigen.all[2,]<1))[1]*4,"-",
                           100- as.numeric(table(eigenall.order < 2.1))[1]*4,"%")),pos=2)
mtext(side=2,line=3,"Eigenvalues")
mtext(side=1,line=2.5,"Number of Dimensions")
axis(side=1,at=seq(2,10,2),labels=seq(2,10,2))
axis(side=1,at=c(1,3),labels=c(1,3))
abline(h=1,lty=2,col="red")

par(mar=c(4,1,4,4))
plot(eigen.rcv[,1],xlim=c(1,10),pch=16,ylim=c(0,20),type='l',ylab="",xlab="",yaxt='n',col="grey60",lwd=1,xaxt='n')
abline(v=2,col="black",lwd=4)
for(i in 1:ncol(eigen.rcv)){
  lines(eigen.rcv[,i],col="grey60",lwd=1)
}
axis(side=4,at=seq(0,20,5),labels=seq(0,20,5),las=1)
mtext(side=3,line=2.25,"(b) Observed Roll Call Votes",cex=1.5)
mtext(side=3,line=.75,print(paste("Average % RCV  = ",round(mean(sim.info.all$percent.rcv),0),"%")),cex=1)
text(x=10,y=19,print(paste("k = ",k)),cex=1,pos=2)
text(x=10,y=15,print(paste0("Missed Second\n Dimension: ",
                            100-as.numeric(table(eigen.rcv[2,]<1))[1]*4,"-",
                            100-as.numeric(table(eigenrcv.order < 2.1))[1]*4,"%")),pos=2)
mtext(side=4,line=3,"Eigenvalues")
mtext(side=1,line=2.5,"Number of Dimensions")
axis(side=1,at=seq(2,10,2),labels=seq(2,10,2))
axis(side=1,at=c(1,3),labels=c(1,3))
abline(h=1,lty=2,col="red")
#dev.off()

#column 4
#load("table-oae3-data-col4.RData")
mean(sim.info.all$total)
mean(sim.info.all$rcv)/mean(sim.info.all$total)
mean(sim.info.all2$percent.all.d1,na.rm=TRUE)
mean(sim.info.all2$percent.all.d2,na.rm=TRUE)
as.numeric(table(eigenall.order < 2.1))[1]*4
as.numeric(table(eigen.all[2,]<1))[1]*4 
mean(sim.info.all2$percent.rcv.d1,na.rm=TRUE)
mean(sim.info.all2$percent.rcv.d2,na.rm=TRUE)
as.numeric(table(eigenrcv.order < 2.1))[1]*4
as.numeric(table(eigen.rcv[2,]<1))[1]*4 

#Generate Figure OA.E13, Row 3
#pdf(".pdf",width=8,height=4)
par(mfrow=c(1,2),mar=c(4,4,4,1))
plot(eigen.all[,1],xlim=c(1,10),ylim=c(0,20),type='l',ylab="",xlab="",yaxt='n',col="grey60",lwd=1,xaxt='n')
abline(v=2,col="black",lwd=4)
for(i in 1:ncol(eigen.all)){
  lines(eigen.all[,i],col="grey60",lwd=1)
}
axis(side=2,at=seq(0,20,5),labels=seq(0,20,5),las=1)
mtext(side=3,line=2.25,"(a) All Votes",cex=1.5)
mtext(side=3,line=.75,print(paste("Average Total Votes  = ",round(mean(sim.info.all$total),0))),cex=1)
text(x=10,y=19,print(paste("k = ",k)),cex=1,pos=2)
text(x=10,y=15,print(paste("Missed Second\n Dimension: ",
                           100-as.numeric(table(eigen.all[2,]<1))[1]*4,"-",
                           100- as.numeric(table(eigenall.order < 2.85))[1]*4,"%")),pos=2)
mtext(side=2,line=3,"Eigenvalues")
mtext(side=1,line=2.5,"Number of Dimensions")
axis(side=1,at=seq(2,10,2),labels=seq(2,10,2))
axis(side=1,at=c(1,3),labels=c(1,3))
abline(h=1,lty=2,col="red")

par(mar=c(4,1,4,4))
plot(eigen.rcv[,1],xlim=c(1,10),pch=16,ylim=c(0,20),type='l',ylab="",xlab="",yaxt='n',col="grey60",lwd=1,xaxt='n')
abline(v=2,col="black",lwd=4)
for(i in 1:ncol(eigen.rcv)){
  lines(eigen.rcv[,i],col="grey60",lwd=1)
}
axis(side=4,at=seq(0,20,5),labels=seq(0,20,5),las=1)
mtext(side=3,line=2.25,"(b) Observed Roll Call Votes",cex=1.5)
mtext(side=3,line=.75,print(paste("Average % RCV  = ",round(mean(sim.info.all$percent.rcv),0),"%")),cex=1)
text(x=10,y=19,print(paste("k = ",k)),cex=1,pos=2)
text(x=10,y=15,print(paste0("Missed Second\n Dimension: ",
                            100-as.numeric(table(eigen.rcv[2,]<1))[1]*4,"-",
                            100-as.numeric(table(eigenrcv.order < 2.85))[1]*4,"%")),pos=2)
mtext(side=4,line=3,"Eigenvalues")
mtext(side=1,line=2.5,"Number of Dimensions")
axis(side=1,at=seq(2,10,2),labels=seq(2,10,2))
axis(side=1,at=c(1,3),labels=c(1,3))
abline(h=1,lty=2,col="red")
#dev.off()


#column 5
#load("table-oae3-data-col5.RData")
mean(sim.info.all$total)
mean(sim.info.all$rcv)/mean(sim.info.all$total)
mean(sim.info.all2$percent.all.d1,na.rm=TRUE)
mean(sim.info.all2$percent.all.d2,na.rm=TRUE)
as.numeric(table(eigenall.order < 2.1))[1]*4
as.numeric(table(eigen.all[2,]<1))[1]*4 
mean(sim.info.all2$percent.rcv.d1,na.rm=TRUE)
mean(sim.info.all2$percent.rcv.d2,na.rm=TRUE)
as.numeric(table(eigenrcv.order < 2.1))[1]*4
as.numeric(table(eigen.rcv[2,]<1))[1]*4 

#Generate Figure OA.E14, Row 1
#pdf(".pdf",width=8,height=4)
par(mfrow=c(1,2),mar=c(4,4,4,1))
plot(eigen.all[,1],xlim=c(1,10),ylim=c(0,20),type='l',ylab="",xlab="",yaxt='n',col="grey60",lwd=1,xaxt='n')
abline(v=2,col="black",lwd=4)
for(i in 1:ncol(eigen.all)){
  lines(eigen.all[,i],col="grey60",lwd=1)
}
axis(side=2,at=seq(0,20,5),labels=seq(0,20,5),las=1)
mtext(side=3,line=2.25,"(a) All Votes",cex=1.5)
mtext(side=3,line=.75,print(paste("Average Total Votes  = ",round(mean(sim.info.all$total),0))),cex=1)
text(x=10,y=19,print(paste("k = ",k)),cex=1,pos=2)
text(x=10,y=15,print(paste("Missed Second\n Dimension: ",
                           100- as.numeric(table(eigenall.order < 2))[1]*4,"%")),pos=2)
mtext(side=2,line=3,"Eigenvalues")
mtext(side=1,line=2.5,"Number of Dimensions")
axis(side=1,at=seq(2,10,2),labels=seq(2,10,2))
axis(side=1,at=c(1,3),labels=c(1,3))
abline(h=1,lty=2,col="red")

par(mar=c(4,1,4,4))
plot(eigen.rcv[,1],xlim=c(1,10),pch=16,ylim=c(0,20),type='l',ylab="",xlab="",yaxt='n',col="grey60",lwd=1,xaxt='n')
abline(v=2,col="black",lwd=4)
for(i in 1:ncol(eigen.rcv)){
  lines(eigen.rcv[,i],col="grey60",lwd=1)
}
axis(side=4,at=seq(0,20,5),labels=seq(0,20,5),las=1)
mtext(side=3,line=2.25,"(b) Observed Roll Call Votes",cex=1.5)
mtext(side=3,line=.75,print(paste("Average % RCV  = ",round(mean(sim.info.all$percent.rcv),0),"%")),cex=1)
text(x=10,y=19,print(paste("k = ",k)),cex=1,pos=2)
text(x=10,y=15,print(paste0("Missed Second\n Dimension: ",
                            100-as.numeric(table(eigen.rcv[2,]<1))[1]*4,"-",
                            100-as.numeric(table(eigenrcv.order < 2))[1]*4,"%")),pos=2)
mtext(side=4,line=3,"Eigenvalues")
mtext(side=1,line=2.5,"Number of Dimensions")
axis(side=1,at=seq(2,10,2),labels=seq(2,10,2))
axis(side=1,at=c(1,3),labels=c(1,3))
abline(h=1,lty=2,col="red")
#dev.off()


#column 6
#load("table-oae3-data-col6.RData")
mean(sim.info.all$total)
mean(sim.info.all$rcv)/mean(sim.info.all$total)
mean(sim.info.all2$percent.all.d1,na.rm=TRUE)
mean(sim.info.all2$percent.all.d2,na.rm=TRUE)
as.numeric(table(eigenall.order < 2.1))[1]*4
as.numeric(table(eigen.all[2,]<1))[1]*4 
mean(sim.info.all2$percent.rcv.d1,na.rm=TRUE)
mean(sim.info.all2$percent.rcv.d2,na.rm=TRUE)
as.numeric(table(eigenrcv.order < 2.1))[1]*4
as.numeric(table(eigen.rcv[2,]<1))[1]*4 

#Generate Figure OA.E14, Row 2
#pdf(".pdf",width=8,height=4)
par(mfrow=c(1,2),mar=c(4,4,4,1))
plot(eigen.all[,1],xlim=c(1,10),ylim=c(0,20),type='l',ylab="",xlab="",yaxt='n',col="grey60",lwd=1,xaxt='n')
abline(v=2,col="black",lwd=4)
for(i in 1:ncol(eigen.all)){
  lines(eigen.all[,i],col="grey60",lwd=1)
}
axis(side=2,at=seq(0,20,5),labels=seq(0,20,5),las=1)
mtext(side=3,line=2.25,"(a) All Votes",cex=1.5)
mtext(side=3,line=.75,print(paste("Average Total Votes  = ",round(mean(sim.info.all$total),0))),cex=1)
text(x=10,y=19,print(paste("k = ",k)),cex=1,pos=2)
text(x=10,y=15,print(paste("Missed Second\n Dimension: ",
                           100- as.numeric(table(eigenall.order < 2))[1]*4,"%")),pos=2)
mtext(side=2,line=3,"Eigenvalues")
mtext(side=1,line=2.5,"Number of Dimensions")
axis(side=1,at=seq(2,10,2),labels=seq(2,10,2))
axis(side=1,at=c(1,3),labels=c(1,3))
abline(h=1,lty=2,col="red")

par(mar=c(4,1,4,4))
plot(eigen.rcv[,1],xlim=c(1,10),pch=16,ylim=c(0,20),type='l',ylab="",xlab="",yaxt='n',col="grey60",lwd=1,xaxt='n')
abline(v=2,col="black",lwd=4)
for(i in 1:ncol(eigen.rcv)){
  lines(eigen.rcv[,i],col="grey60",lwd=1)
}
axis(side=4,at=seq(0,20,5),labels=seq(0,20,5),las=1)
mtext(side=3,line=2.25,"(b) Observed Roll Call Votes",cex=1.5)
mtext(side=3,line=.75,print(paste("Average % RCV  = ",round(mean(sim.info.all$percent.rcv),0),"%")),cex=1)
text(x=10,y=19,print(paste("k = ",k)),cex=1,pos=2)
text(x=10,y=15,print(paste0("Missed Second\n Dimension: ",
                            100-as.numeric(table(eigen.rcv[2,]<1))[1]*4,"-",
                            100-as.numeric(table(eigenrcv.order < 2))[1]*4,"%")),pos=2)
mtext(side=4,line=3,"Eigenvalues")
mtext(side=1,line=2.5,"Number of Dimensions")
axis(side=1,at=seq(2,10,2),labels=seq(2,10,2))
axis(side=1,at=c(1,3),labels=c(1,3))
abline(h=1,lty=2,col="red")
#dev.off()


###########################################################
##                     Figure OA.E15                     ##
###########################################################

#load("fig-oae15-data.RData")

#Generate Figure OA.E15
#pdf(".pdf",width=8,height=4)
par(mfrow=c(1,2),mar=c(4,4,4,1))
plot(eigen.all[,1],xlim=c(1,10),ylim=c(0,20),type='l',ylab="",xlab="",yaxt='n',col="grey60",lwd=1,xaxt='n')
abline(v=2,col="black",lwd=4)
for(i in 1:ncol(eigen.all)){
  lines(eigen.all[,i],col="grey60",lwd=1)
}
axis(side=2,at=seq(0,20,5),labels=seq(0,20,5),las=1)
mtext(side=3,line=2.25,"(a) All Votes",cex=1.5)
mtext(side=3,line=.75,print(paste("Average Total Votes  = ",round(mean(sim.info.all$total),0))),cex=1)
text(x=10,y=19,print(paste("k = ",k)),cex=1,pos=2)
text(x=10,y=15,print(paste("Missed Second\n Dimension: ",
                           100-as.numeric(table(eigen.all[2,]<1))[1]*4,"-",
                           100- as.numeric(table(eigenall.order < 2.1))[1]*4,"%")),pos=2)
mtext(side=2,line=3,"Eigenvalues")
mtext(side=1,line=2.5,"Number of Dimensions")
axis(side=1,at=seq(2,10,2),labels=seq(2,10,2))
axis(side=1,at=c(1,3),labels=c(1,3))
abline(h=1,lty=2,col="red")

par(mar=c(4,1,4,4))
plot(eigen.rcv[,1],xlim=c(1,10),pch=16,ylim=c(0,20),type='l',ylab="",xlab="",yaxt='n',col="grey60",lwd=1,xaxt='n')
abline(v=2,col="black",lwd=4)
for(i in 1:ncol(eigen.rcv)){
  lines(eigen.rcv[,i],col="grey60",lwd=1)
}
axis(side=4,at=seq(0,20,5),labels=seq(0,20,5),las=1)
mtext(side=3,line=2.25,"(b) Observed Roll Call Votes",cex=1.5)
mtext(side=3,line=.75,print(paste("Average % RCV  = ",round(mean(sim.info.all$percent.rcv),0),"%")),cex=1)
text(x=10,y=19,print(paste("k = ",k)),cex=1,pos=2)
text(x=10,y=15,print(paste0("Missed Second\n Dimension: ",
                            100-as.numeric(table(eigen.rcv[2,]<1))[1]*4,"-",
                            100-as.numeric(table(eigenrcv.order < 2.1))[1]*4,"%")),pos=2)
mtext(side=4,line=3,"Eigenvalues")
mtext(side=1,line=2.5,"Number of Dimensions")
axis(side=1,at=seq(2,10,2),labels=seq(2,10,2))
axis(side=1,at=c(1,3),labels=c(1,3))
abline(h=1,lty=2,col="red")
#dev.off()


###########################################################
##                      Table OA.E4                      ##
###########################################################

#column 1
#load("table-oae4-data-col1.RData")
mean(sim.info.all$total)
mean(sim.info.all$rcv / (sim.info.all$non.rcv+sim.info.all$rcv)) 
mean(sim.info.all2$percent.all.d1,na.rm=TRUE)
mean(sim.info.all2$percent.all.d2,na.rm=TRUE)
100-(25-as.numeric(table(eigen.all[2,]<2.1))[1])*4
100-(25-as.numeric(table(eigen.all[2,]<1))[1])*4
mean(sim.info.all2$percent.rcv.d1,na.rm=TRUE)
mean(sim.info.all2$percent.rcv.d2,na.rm=TRUE)
100-(25-as.numeric(table(eigen.rcv[2,]<2.1))[1])*4
100-(25-as.numeric(table(eigen.rcv[2,]<1))[1])*4

#column 2
#load("table-oae4-data-col2.RData")
mean(sim.info.all$total)
mean(sim.info.all$rcv / (sim.info.all$non.rcv+sim.info.all$rcv)) 
mean(sim.info.all2$percent.all.d1,na.rm=TRUE)
mean(sim.info.all2$percent.all.d2,na.rm=TRUE)
100-(25-as.numeric(table(eigen.all[2,]<2.1))[1])*4
100-(25-as.numeric(table(eigen.all[2,]<1))[1])*4
mean(sim.info.all2$percent.rcv.d1,na.rm=TRUE)
mean(sim.info.all2$percent.rcv.d2,na.rm=TRUE)
100-(25-as.numeric(table(eigen.rcv[2,]<2.1))[1])*4
100-(25-as.numeric(table(eigen.rcv[2,]<1))[1])*4

#column 3
#load("table-oae4-data-col3.RData")
mean(sim.info.all$total)
mean(sim.info.all$rcv / (sim.info.all$non.rcv+sim.info.all$rcv)) 
mean(sim.info.all2$percent.all.d1,na.rm=TRUE)
mean(sim.info.all2$percent.all.d2,na.rm=TRUE)
100-(25-as.numeric(table(eigen.all[2,]<2.1))[1])*4
100-(25-as.numeric(table(eigen.all[2,]<1))[1])*4
mean(sim.info.all2$percent.rcv.d1,na.rm=TRUE)
mean(sim.info.all2$percent.rcv.d2,na.rm=TRUE)
100-(25-as.numeric(table(eigen.rcv[2,]<2.1))[1])*4
100-(25-as.numeric(table(eigen.rcv[2,]<1))[1])*4

#column 4
#load("table-oae4-data-col4.RData")
mean(sim.info.all$total)
mean(sim.info.all$rcv / (sim.info.all$non.rcv+sim.info.all$rcv)) 
mean(sim.info.all2$percent.all.d1,na.rm=TRUE)
mean(sim.info.all2$percent.all.d2,na.rm=TRUE)
100-(25-as.numeric(table(eigen.all[2,]<2.1))[1])*4
100-(25-as.numeric(table(eigen.all[2,]<1))[1])*4
mean(sim.info.all2$percent.rcv.d1,na.rm=TRUE)
mean(sim.info.all2$percent.rcv.d2,na.rm=TRUE)
100-(25-as.numeric(table(eigen.rcv[2,]<2.1))[1])*4
100-(25-as.numeric(table(eigen.rcv[2,]<1))[1])*4


###########################################################
##                     Figure OA.E16                     ##
###########################################################

library(wnominate)
library(pscl)

alpha <- .5
k <- .35
N <- 100
v <- 500
nsim <- 1

eigen.all <- NULL
eigen.rcv <- NULL
sim.info.all <- NULL

for(j in 1:nsim){
  phi <- runif(1,.2,.8)
  b1 <- runif(v,0,1)
  b2 <- runif(v,0,1)
  sq1 <- runif(v,0,1)
  sq2 <- runif(v,0,1)
  
  npL <- round((1-phi)*N,digits=0)
  npR <- round((phi)*N,digits=0)
  d1 <- runif(1,.5,2)
  d2 <- runif(1,.5,2)
  
    hypo.pro.1d <- cbind(rep(1,length(b1)),rep(d1,length(b1)),b1,sq1)
    colnames(hypo.pro.1d) <- c("dim","d","b","sq")
    hypo.pro.2d <- cbind(rep(2,length(b2)),rep(d2,length(b2)),b2,sq2)
    colnames(hypo.pro.2d) <- c("dim","d","b","sq")
    hypo.pro.all <- as.data.frame(rbind(hypo.pro.1d,hypo.pro.2d))
    hypo.pro.all$id <- seq(1,nrow(hypo.pro.all),1)
    
    #ideal points
    xL1 <- c(0,runif(npL-1,min=0-d1,max=d1)) #Going to use this that i've put the party leaders in the first and last spots
    xL2 <- c(0,runif(npL-1,min=0-d2,max=d2))
    xR1 <- c(runif(npR-1,min=1-d1,max=1+d1),1) #Going to use this that i've put the party leaders in the first and last spots
    xR2 <- c(runif(npR-1,min=1-d2,max=1+d2),1)
    left.party <- as.data.frame(cbind(xL1,xL2,rep(0,npL)))
    colnames(left.party) <- c("x1","x2","partyid")
    right.party <- as.data.frame(cbind(xR1,xR2,rep(1,npR)))
    colnames(right.party) <- c("x1","x2","partyid")
    all.leg <- as.data.frame(rbind(left.party,right.party))
    
    #who votes for it?
    vote.func1 <- function(x1,b1,sq1){
      as.numeric(abs(x1-b1)  < abs(x1-sq1) )
    }
    vote.func2 <- function(x2,b2,sq2){
      as.numeric(abs(x2-b2)  < abs(x2-sq2) )
    }
    out <- matrix(data=NA,nrow=N,ncol=v*2);dim(out)
    for(i in 1:N){
      for(j in 1:v){
        out[i,j] <- vote.func1(all.leg$x1[i],b1[j],sq1[j])
      }
    }
    
    for(i in 1:N){
      for(j in (v+1):(2*v)){
        out[i,j] <- vote.func2(all.leg$x2[i],b2[j-v],sq2[j-v])
      }
    }
    
    #pass?    
    passes <- as.numeric(colSums(out) > .5*N)
    hypo.pro.all <- as.data.frame(cbind(hypo.pro.all,passes))
    
    for (i in 1:nrow(hypo.pro.all)){
      hypo.pro.all$outcome[i] <- if(hypo.pro.all$passes[i]==0) hypo.pro.all$sq[i] else hypo.pro.all$b[i]
    }
    
    left.aligned <- rep(NA,length=ncol(out))
    for(i in 1:ncol(out)){
      left.aligned[i] <- round((sum(as.numeric(out[2:npL,i]==out[1,i])))/npL,4)
    }
    right.aligned <- rep(NA,length=ncol(out))
    for(i in 1:ncol(out)){
      right.aligned[i] <- round((sum(as.numeric(out[(npL+1):(N-1),i]==out[N,i])))/npR,4)
    }
    
    out.cohesion <- matrix(data=NA,nrow=N,ncol=v*2)
    for(j in 1:(v*2)){
      for(i in 1:npL){
        out.cohesion[i,j] <- if (out[i,j] == out[1,j]) {
          left.aligned[j] - right.aligned[j]
        } else {
          0
        }
      }
    }
    for(j in 1:(v*2)){
      for(i in (npL+1):N){
        out.cohesion[i,j] <- if (out[i,j] == out[N,j]) {
          right.aligned[j]-left.aligned[j]
        } else {
          0
        }
      }
    }
    out.cohesion.abs <- matrix(data=NA,nrow=N,ncol=v*2)
    for(j in 1:(v*2)){
      for(i in 1:npL){
        out.cohesion.abs[i,j] <- if (out[i,j] == out[1,j]) {
          left.aligned[j] 
        } else {
          0
        }
      }
    }
    for(j in 1:(v*2)){
      for(i in (npL+1):N){
        out.cohesion.abs[i,j] <- if (out[i,j] == out[N,j]) {
          right.aligned[j]
        } else {
          0
        }
      }
    }
    
    #roll call votes requested?
    out.cohesion.kmat <- matrix(data=NA,nrow=N,ncol=v*2)
    for(i in 1:N){
      for(j in 1:(v*2)){
        out.cohesion.kmat[i,j] <- as.numeric(out.cohesion[i,j] > k)
      }
    }
    hypo.pro.all$left.cscore <- out.cohesion[1,] 
    hypo.pro.all$right.cscore <- out.cohesion[N,]
    hypo.pro.all$left.cscore.abs <- out.cohesion.abs[1,] 
    hypo.pro.all$right.cscore.abs <- out.cohesion.abs[N,]
    hypo.pro.all$left.rcv.request <- out.cohesion.kmat[1,]
    hypo.pro.all$right.rcv.request <- out.cohesion.kmat[N,]    
    hypo.pro.all$rcv.request <- as.numeric(colSums(out.cohesion.kmat)>0)
    
    all.proposals <- hypo.pro.all
    rcv.proposals <- hypo.pro.all[which(hypo.pro.all$rcv.request==1),]
    vote.matrix.all <- out[,all.proposals$id]
    vote.matrix.rcv <- out[,rcv.proposals$id]
    
    #Dimensionality
    votemat.all <- rollcall(vote.matrix.all, yea=1, nay=0)
    votemat.rcv <- rollcall(vote.matrix.rcv, yea=1, nay=0)
    
    res.all <- wnominate(votemat.all,polarity = c(100,100),dims=2,
                         verbose=FALSE)
    res.rcv <- wnominate(votemat.rcv,polarity = c(100,100),dims=2,
                         verbose=FALSE)
    eigen.all <- cbind(eigen.all,res.all$eigenvalues)
    eigen.rcv <- cbind(eigen.rcv,res.rcv$eigenvalues)
    
    temp.all <- as.numeric(table(all.proposals$d))
    temp.rcv <- as.numeric(table(rcv.proposals$d))
    sim.info <- c(nrow(all.proposals),nrow(rcv.proposals),phi,d1,d2,temp.all[2],temp.all[1],temp.rcv[2],temp.rcv[1])
    sim.info.all <- rbind(sim.info.all,sim.info)
}

eigen.all <- eigen.all[,1:25]
eigen.rcv <- eigen.rcv[,1:25]
colnames(sim.info.all) <- c("non.rcv","rcv","phi","d1","d2","all.votes.d1","all.votes.d2","rcv.votes.d1","rcv.votes.d2")
sim.info.all <- as.data.frame(sim.info.all)
sim.info.all <- sim.info.all[1:25,]
sim.info.all$total <- sim.info.all$non.rcv+sim.info.all$rcv
sim.info.all$percent.rcv <- (sim.info.all$rcv / sim.info.all$total)*100

#load("fig-oae16-data.RData")
#Generate Figure OA.E16
#pdf(".pdf",width=8,height=4)
par(mfrow=c(1,2),mar=c(4,4,4,1))
plot(eigen.all[,1],xlim=c(1,10),ylim=c(0,20),type='l',ylab="",xlab="",yaxt='n',col="grey60",lwd=1,xaxt='n')
abline(v=2,col="black",lwd=4)
for(i in 1:ncol(eigen.all)){
  lines(eigen.all[,i],col="grey60",lwd=1)
}
axis(side=2,at=seq(0,20,5),labels=seq(0,20,5),las=1)
mtext(side=3,line=2.25,"(a) All Votes",cex=1.5)
mtext(side=3,line=.75,print(paste("Average Total Votes  = 1000")),cex=1)
text(x=10,y=19,print(paste("k = ",k)),cex=1,pos=2)
text(x=10,y=15,print(paste("Missed Second\n Dimension: ","0-8%" )),pos=2)
mtext(side=2,line=3,"Eigenvalues")
mtext(side=1,line=2.5,"Number of Dimensions")
axis(side=1,at=seq(2,10,2),labels=seq(2,10,2))
axis(side=1,at=c(1,3),labels=c(1,3))
abline(h=1,lty=2,col="red")

par(mar=c(4,1,4,4))
plot(eigen.rcv[,1],xlim=c(1,10),pch=16,ylim=c(0,20),type='l',ylab="",xlab="",yaxt='n',col="grey60",lwd=1,xaxt='n')
abline(v=2,col="black",lwd=4)
for(i in 1:ncol(eigen.rcv)){
  lines(eigen.rcv[,i],col="grey60",lwd=1)
}
axis(side=4,at=seq(0,20,5),labels=seq(0,20,5),las=1)
mtext(side=3,line=2.25,"(b) Observed Roll Call Votes",cex=1.5)
mtext(side=3,line=.75,print(paste("Average % RCV  = ",round(mean(sim.info.all$percent.rcv),0),"%")),cex=1)
text(x=10,y=19,print(paste("k = ",k)),cex=1,pos=2)
text(x=10,y=15,print(paste0("Missed Second\n Dimension: ",((as.numeric(table(eigen.rcv[2,]<=1)[2])/25)*100),"-40%" )),pos=2)
mtext(side=4,line=3,"Eigenvalues")
mtext(side=1,line=2.5,"Number of Dimensions")
axis(side=1,at=seq(2,10,2),labels=seq(2,10,2))
axis(side=1,at=c(1,3),labels=c(1,3))
abline(h=1,lty=2,col="red")
#dev.off()